Statistical analysis

Litter decomposition

Author

Marcelo Araya-Salas & Andrea Vincent

Published

February 17, 2023

Source code and data found at

Purpose

Evaluate role of nutrient availability on litter decomposition

1 Nutrient change

1.1 Nitrogen

Content

Code

nutr <- read.csv("./data/raw/litter-nutrients-mas.csv")

nutr$plot.f <- as.factor(nutr$plot)

nutr$days.sc <- scale(nutr$days)

nutr$litter.n.content.prop.initial <- nutr$litter.n.content.perc.initial/100

# remove plot 9
sub.nutr <- nutr[nutr$plot != 9, ]

agg_n <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, mean)

agg_n$sd <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, sd)$litter.n.content.perc.initial

agg_n$se <- aggregate(litter.n.content.perc.initial ~ colecta + treat +
    days, nutr, se)$litter.n.content.perc.initial

agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_n, aes(x = days, y = litter.n.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.n.content.perc.initial +
    se, ymin = litter.n.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter N content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_n$days), labels = unique(agg_n$days)) +
    theme(legend.position = c(0.9, 0.8))

Code

fit.n <- brm(litter.n.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, cores = chains, family = Beta(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    prior = prior, file = "./data/processed/litter.n.content.rem_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/litter.n.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.2 litter.n.content.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	litter.n.content.prop.initial ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	19631.06	22594.55	272268972

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.340	-0.634	-0.048	1	27796.92	26514.82
treatNP	-0.051	-0.344	0.238	1	27673.41	24838.46
treatP	-0.060	-0.349	0.231	1	28095.99	25183.40
days.sc	-0.503	-0.712	-0.307	1	19631.06	22594.55
treatN:days.sc	-0.094	-0.369	0.183	1	25338.01	27759.09
treatNP:days.sc	0.039	-0.227	0.314	1	25217.60	27081.45
treatP:days.sc	0.048	-0.220	0.323	1	23691.70	28287.99

Concentration

Code

agg_n <- aggregate(perc.n ~ colecta + treat + days, nutr, mean)

agg_n$sd <- aggregate(perc.n ~ colecta + treat + days, nutr, sd)$perc.n

agg_n$se <- aggregate(perc.n ~ colecta + treat + days, nutr, se)$perc.n

agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_n, aes(x = days, y = perc.n, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.n + se, ymin = perc.n -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_n$days),
    labels = unique(agg_n$days)) + theme(legend.position = c(0.3,
    0.8))

Code

nutr$prop.n <- nutr$perc.n/100

fit.perc.n <- brm(prop.n ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, core = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.n.perc_model", file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/litter.n.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.3 litter.n.perc_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.n ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	16686.78	23193.79	757288462

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.097	-0.240	0.046	1	25471.79	25380.38
treatNP	-0.052	-0.192	0.087	1	25656.85	26419.36
treatP	-0.025	-0.163	0.114	1	24461.38	24258.37
days.sc	0.134	0.041	0.227	1	16686.78	23193.79
treatN:days.sc	-0.075	-0.205	0.052	1	20777.70	25856.35
treatNP:days.sc	-0.030	-0.160	0.099	1	20665.11	26156.06
treatP:days.sc	-0.012	-0.138	0.116	1	20899.42	26866.06

1.4 Phosphorus

Content

Code

nutr$litter.p.content.prop.initial <- nutr$litter.p.content.perc.initial/100

# excluding plot 9
agg_p <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, mean)

agg_p$sd <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, sd)$litter.p.content.perc.initial

agg_p$se <- aggregate(litter.p.content.perc.initial ~ colecta + treat +
    days, sub.nutr, se)$litter.p.content.perc.initial


agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_p, aes(x = days, y = litter.p.content.perc.initial, color = treat)) +
    geom_point(size = 2, position = pd) + geom_errorbar(aes(ymax = litter.p.content.perc.initial +
    se, ymin = litter.p.content.perc.initial - se), width = 0, position = pd) +
    geom_line(size = 1.2, position = pd) + scale_color_viridis_d(alpha = 0.5) +
    labs(x = "Time (days)", y = "Litter P content (% initial)", color = "Treatment") +
    scale_x_continuous(breaks = unique(agg_p$days), labels = unique(agg_p$days)) +
    theme(legend.position = c(0.9, 0.8))

Code

fit.p <- brm(litter.p.content.prop.initial ~ treat * days.sc + (1 |
    plot.f), data = nutr, chains = chains, prior = prior, cores = chains,
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/litter.p.content.rem_model", file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/litter.p.content.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.5 litter.p.content.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	litter.p.content.prop.initial ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	13618.61	17352.52	2063397055

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.090	-1.258	1.081	1	13662.76	18828.78
treatNP	0.009	-1.148	1.172	1	13618.61	17352.52
treatP	0.754	-0.399	1.907	1	14549.30	18770.67
days.sc	-0.134	-0.473	0.205	1	17312.08	23553.03
treatN:days.sc	0.001	-0.459	0.459	1	20983.83	27133.63
treatNP:days.sc	0.021	-0.441	0.482	1	21015.75	25715.30
treatP:days.sc	-0.168	-0.620	0.288	1	20970.54	26482.58

Concentration

Code

agg_p <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, mean)

agg_p$sd <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, sd)$perc.p

agg_p$se <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, se)$perc.p

agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))

pd <- position_dodge(15)

ggplot(agg_p, aes(x = days, y = perc.p, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = perc.p + se, ymin = perc.p -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_p$days),
    labels = unique(agg_p$days)) + theme(legend.position = c(0.9,
    0.8))

Code

# convert to proportions to use beta distribution
sub.nutr$prop.p <- sub.nutr$perc.p/100

fit.perc.p <- brm(prop.p ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, cores = chains, family = Beta(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), prior = prior,
    file = "./data/processed/litter.p.perc_model", file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/litter.p.perc_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.6 litter.p.perc_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.p ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	17993.25	23361.17	1797840914

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	-0.106	-0.634	0.409	1	28569.54	28049.17
treatNP	0.031	-0.475	0.530	1	27422.98	28348.78
treatP	0.077	-0.468	0.610	1	28087.50	28163.62
days.sc	-0.034	-0.388	0.302	1	17993.25	23361.17
treatN:days.sc	-0.001	-0.486	0.478	1	22840.20	28837.19
treatNP:days.sc	0.040	-0.430	0.503	1	21964.87	27252.32
treatP:days.sc	0.028	-0.475	0.531	1	23659.98	28583.96

Takeaways

Litter P content is significantly lower in plus N treatment plots than in control plot, after accounting for variation explained by time

1.7 Remaining litter

Code

fit <- brm(prop.litter.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model",
    file_refit = "on_change")

dat$days.fc <- as.numeric(as.factor(dat$days))

# monotonic effect of time
fit_mo <- brm(prop.litter.rem ~ trat * mo(days.fc) + (1 | plot.f),
    data = dat, chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model_monotonic",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/prop.litter.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.8 prop.litter.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.litter.rem ~ trat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	17979.85	23910.9	1631599382

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
tratN	-0.038	-0.203	0.125	1	30540.12	29362.46
tratNP	0.074	-0.085	0.233	1	29883.05	29157.97
tratP	0.032	-0.127	0.192	1	30065.88	30156.88
days.sc	-0.683	-0.800	-0.570	1	17979.85	23910.90
tratN:days.sc	-0.047	-0.207	0.113	1	23277.78	30456.01
tratNP:days.sc	0.106	-0.045	0.259	1	22364.84	28848.92
tratP:days.sc	0.056	-0.097	0.212	1	22771.51	27130.32

Code

html_summary(read.file = "./data/processed/prop.litter.rem_model_monotonic.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.9 prop.litter.rem_model_monotonic

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	prop.litter.rem ~ trat * mo(days.fc) + (1 \| plot.f)	20000	4	1	10000	0	0	17641.98	23834.62	407718541

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
tratN	0.038	-0.189	0.270	1	19420.68	24276.66
tratNP	-0.037	-0.259	0.187	1	17828.65	24753.65
tratP	-0.009	-0.238	0.221	1	17641.98	23834.62

1.10 Remaining wood

Code

dat$prop.wood.rem <- dat$perc.wood.rem/100

fit2 <- brm(prop.wood.rem ~ trat * days.sc + (1 | plot.f), data = dat,
    chains = chains, family = Beta(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/prop.wood.rem_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/prop.wood.rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

1.11 prop.wood.rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.wood.rem ~ trat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	16872.71	23764.8	1828689869

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
tratN	0.118	-0.128	0.366	1	25420.93	24370.70
tratNP	0.254	0.012	0.496	1	24798.43	26473.85
tratP	0.197	-0.045	0.438	1	25725.64	25244.62
days.sc	-0.910	-1.094	-0.737	1	16872.71	23764.80
tratN:days.sc	-0.006	-0.244	0.234	1	20661.18	27810.64
tratNP:days.sc	0.097	-0.140	0.333	1	20685.17	28019.73
tratP:days.sc	0.175	-0.058	0.412	1	20811.62	27760.94

2 K

2.1 Litter

Correlation between Silvia’s and Andrea’s K

Code

k_vals <- read.csv("./data/raw/k-values-corr.csv")

cor(k_vals$av.k.litter, k_vals$sil.k.litter)

[1] 0.9327804

Code

cor(k_vals$av.k.wood, k_vals$sil.k.wood)

[1] 0.9621473

Comparing all treatments vs control

Code

fit_k.litter <- brm(sil.k.litter ~ Treatment + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_model",
    file_refit = "on_change")

Code

html_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "",
    highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_model.rds")

2.2 k_litter_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ Treatment + (1 \| quadrat)	20000	4	1	10000	5	0	22486.82	24586.15	224004760

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
TreatmentN	0.238	-0.176	0.659	1	22486.82	24586.15
TreatmentNP	-0.156	-0.572	0.261	1	25034.00	26992.34
TreatmentP	-0.006	-0.425	0.412	1	23833.14	25611.23

Phosphorus vs no-phosphorus

Code

fit_k.litter.p.np <- brm(sil.k.litter ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, , control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_p_nop_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/k_litter_p_nop_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.3 k_litter_p_nop_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ p.treat + (1 \| quadrat)	20000	4	1	10000	8	0	38432.94	27572.89	1236972411

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
p.treatwith.p	-0.201	-0.492	0.092	1	38432.94	27572.89

Nitrogen vs no-nitrogen

Code

fit_k.litter.n.nn <- brm(sil.k.litter ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_litter_n_no_n_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/k_litter_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.4 k_litter_n_no_n_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ n.treat + (1 \| quadrat)	20000	4	1	10000	17	0	37023.21	27505.83	956380015

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
n.treatwith.n	0.045	-0.254	0.345	1	37023.21	27505.83

2.5 Wood

Comparing all treatments vs control

Code

agg_kvals <- aggregate(sil.k.wood ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ Treatment, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ Treatment, k_vals, sd)$sil.k.wood

agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N",
    "P", "NP"))

ggplot(k_vals, aes(x = Treatment, y = sil.k.wood)) + geom_violin(fill = "#21908C4D") +
    labs(y = "Decomposition constant (k)") + geom_pointrange(data = agg_kvals,
    aes(ymin = sil.k.wood - se, ymax = sil.k.wood + se))

Code

fit_k.wood.p.treat <- brm(sil.k.wood ~ Treatment + (1 | quadrat),
    data = k_vals, chains = chains, family = gaussian(), iter = iters,
    control = list(adapt_delta = 0.99, max_treedepth = 15), cores = chains,
    prior = prior, file = "./data/processed/k_wood_treatment_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/k_wood_treatment_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

2.6 k_wood_treatment_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ Treatment + (1 \| quadrat)	20000	4	1	10000	465	0	7013.536	25155.17	1147311351

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
TreatmentN	0.031	-0.537	0.604	1	7013.536	26593.78
TreatmentNP	-0.491	-1.067	0.073	1.001	19096.320	25155.17
TreatmentP	-0.471	-1.045	0.105	1.001	18153.985	25925.35

3 Phosphorus vs no-phosphorus

Code

agg_kvals <- aggregate(sil.k.wood ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ p.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ p.treat, k_vals, sd)$sil.k.wood


agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P",
    "P"))

k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))

ggplot(k_vals, aes(x = p.treat, y = sil.k.wood)) + geom_violin(fill = "#21908C4D") +
    labs(y = "Decomposition constant (k)", x = "P treatment") + geom_pointrange(data = agg_kvals,
    aes(ymin = sil.k.wood - se, ymax = sil.k.wood + se))

Code

fit_k.wood.p.no.p <- brm(sil.k.wood ~ p.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_p_no_p_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/k_wood_p_no_p_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

3.1 k_wood_p_no_p_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ p.treat + (1 \| quadrat)	20000	4	1	10000	4	0	31952.25	26238.54	326631910

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
p.treatwith.p	-0.496	-0.892	-0.1	1	31952.25	26238.54

4 Nitrogen vs no-nitrogen

Code

agg_kvals <- aggregate(sil.k.wood ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ n.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ n.treat, k_vals, sd)$sil.k.wood


agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N",
    "N"))

k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))

ggplot(k_vals, aes(x = n.treat, y = sil.k.wood)) + geom_violin(fill = "#21908C4D") +
    labs(y = "Decomposition constant (k)", x = "N treatment") + geom_pointrange(data = agg_kvals,
    aes(ymin = sil.k.wood - se, ymax = sil.k.wood + se))

Code

fit_k.wood.n.no.n <- brm(sil.k.wood ~ n.treat + (1 | quadrat), data = k_vals,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/k_wood_n_no_n_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/k_wood_n_no_n_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

4.1 k_wood_n_no_n_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ n.treat + (1 \| quadrat)	20000	4	1	10000	7	0	30908.88	27135.72	1683655054

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
n.treatN	0.008	-0.414	0.426	1	30908.88	27135.72

5 Nutrient content by remaining litter mass

5.1 Nitrogen

Code

nutr$prom.hoja.reman.sc <- scale(nutr$prom.hoja.reman)

ggplot(data = nutr, aes(x = prom.hoja.reman, y = perc.n, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter N concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code

fit_pern_by_rem_leave <- brm(perc.n ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.2 n_per_by_leave_rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.n ~ prom.hoja.reman.sc * treat + (1 \| plot)	20000	4	1	10000	0	0	22712.6	26086.89	202030287

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
prom.hoja.reman.sc	-0.152	-0.265	-0.038	1	22712.60	27352.64
treatN	-0.161	-0.345	0.025	1	29404.33	26086.89
treatNP	-0.044	-0.224	0.138	1	29226.14	26624.22
treatP	-0.037	-0.218	0.145	1	28414.57	26316.85
prom.hoja.reman.sc:treatN	0.060	-0.095	0.212	1	27004.70	29358.59
prom.hoja.reman.sc:treatNP	-0.103	-0.269	0.064	1	28381.31	29425.85
prom.hoja.reman.sc:treatP	-0.013	-0.174	0.150	1	27640.20	30237.33

5.3 Phosphorus

Code

ggplot(data = sub.nutr, aes(x = prom.hoja.reman, y = perc.p, color = treat)) +
    geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter P concentration (%)",
    color = "Treatment") + scale_color_viridis_d(alpha = 0.5) + geom_smooth(method = "lm",
    se = FALSE) + scale_x_reverse()

Code

sub.nutr$prom.hoja.reman.sc <- scale(sub.nutr$prom.hoja.reman)

fit_perp_by_rem_leave <- brm(perc.p ~ prom.hoja.reman.sc * treat +
    (1 | plot), data = sub.nutr, chains = chains, family = gaussian(),
    iter = iters, control = list(adapt_delta = 0.99, max_treedepth = 15),
    cores = chains, prior = prior, file = "./data/processed/p_per_by_leave_rem_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/p_per_by_leave_rem_model.rds",
    gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
    remove.intercepts = TRUE)

5.4 p_per_by_leave_rem_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.p ~ prom.hoja.reman.sc * treat + (1 \| plot)	20000	4	1	10000	0	0	18469.62	19550.38	8818990

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
prom.hoja.reman.sc	0.001	-0.002	0.004	1	21082.94	25940.57
treatN	-0.005	-0.012	0.002	1	18785.39	19550.38
treatNP	0.000	-0.007	0.007	1	18469.62	20718.11
treatP	0.003	-0.004	0.011	1	18490.36	20587.53
prom.hoja.reman.sc:treatN	0.000	-0.004	0.005	1	24965.09	28102.51
prom.hoja.reman.sc:treatNP	-0.002	-0.006	0.003	1	25491.25	27830.25
prom.hoja.reman.sc:treatP	-0.001	-0.005	0.004	1	26931.48	27838.27

Takeaways

Nitrogen percentage, but no Phosphorus percentage, increases along with remaining litter mass

6 Nutrient ratios

6.1 C:N

Code

# excluding plot 9
agg_cn <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, mean)

agg_cn$sd <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, sd)$cn.mol.kg

agg_cn$se <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, se)$cn.mol.kg


agg_cn$treat <- factor(agg_cn$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_cn, aes(x = days, y = cn.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cn.mol.kg + se, ymin = cn.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:N ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cn$days),
    labels = unique(agg_cn$days)) + theme(legend.position = c(0.9,
    0.8))

Code

fit.cn <- brm(cn.mol.kg ~ treat * days.sc + (1 | plot.f), data = nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/cn_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/cn_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.2 cn_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	cn.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	23616.14	27188.36	206225990

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	0.309	-1.373	2.015	1	33345.57	27188.36
treatNP	-0.189	-1.881	1.509	1	31802.40	27824.88
treatP	-0.182	-1.863	1.514	1	30821.46	28785.15
days.sc	-3.101	-4.204	-1.996	1	23616.14	27323.41
treatN:days.sc	0.751	-0.752	2.245	1	29469.69	30884.24
treatNP:days.sc	0.439	-1.062	1.942	1	29505.97	31955.60
treatP:days.sc	0.111	-1.391	1.614	1	29297.05	31633.29

6.3 N:P

Code

# excluding plot 9
agg_np <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_np$sd <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$np.mol.kg

agg_np$se <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$np.mol.kg


agg_np$treat <- factor(agg_np$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_np, aes(x = days, y = np.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = np.mol.kg + se, ymin = np.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter N:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_np$days),
    labels = unique(agg_np$days)) + theme(legend.position = c(0.2,
    0.8))

Code

fit.np <- brm(np.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta = 0.99,
        max_treedepth = 15), cores = chains, prior = prior, file = "./data/processed/np_model",
    file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/np_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.4 np_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3)	np.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	23138.62	23571.24	226626333

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	1.211	-3.076	5.572	1	27662.79	23728.89
treatNP	-2.401	-6.601	1.896	1	26490.49	23571.24
treatP	-3.135	-7.696	1.500	1	27059.33	23906.00
days.sc	4.407	1.755	7.072	1	23138.62	26390.49
treatN:days.sc	-1.765	-5.391	1.885	1	26580.63	27306.41
treatNP:days.sc	-1.780	-5.403	1.901	1	27640.74	28593.42
treatP:days.sc	-1.415	-5.318	2.486	1	28596.79	29590.75

6.5 C:P

Code

# excluding plot 9
agg_cp <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    mean)

agg_cp$sd <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    sd)$cp.mol.kg

agg_cp$se <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr,
    se)$cp.mol.kg

agg_cp$treat <- factor(agg_cp$treat, levels = c("C", "N", "P", "NP"))


pd <- position_dodge(15)

ggplot(agg_cp, aes(x = days, y = cp.mol.kg, color = treat)) + geom_point(size = 2,
    position = pd) + geom_errorbar(aes(ymax = cp.mol.kg + se, ymin = cp.mol.kg -
    se), width = 0, position = pd) + geom_line(size = 1.2, position = pd) +
    scale_color_viridis_d(alpha = 0.5) + labs(x = "Time (days)", y = "Litter C:P ratio",
    color = "Treatment") + scale_x_continuous(breaks = unique(agg_cp$days),
    labels = unique(agg_cp$days)) + theme(legend.position = c(0.9,
    0.8))

Code

fit.cp <- brm(cp.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr,
    chains = chains, family = gaussian(), iter = iters, cores = chains,
    prior = prior, file = "./data/processed/cp_model", file_refit = "on_change")

Code

html_summary(read.file = "./data/processed/cp_model.rds", gsub.pattern = "b_treatment|b_",
    gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)

6.6 cp_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3)	cp.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	0	30235.23	26541.32	1179254390

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
treatN	0.386	-19.206	20.050	1	30235.23	27580.58
treatNP	-0.159	-19.492	19.226	1	31461.97	26541.32
treatP	-0.300	-19.898	19.294	1	32550.30	27705.93
days.sc	-6.177	-21.408	9.060	1	30959.55	27830.18
treatN:days.sc	0.319	-17.923	18.305	1	32235.63	27315.59
treatNP:days.sc	-4.800	-23.038	13.518	1	32491.62	27863.26
treatP:days.sc	-2.648	-20.973	15.782	1	32066.82	27837.57

7 Combined model diagnostics

Code

check_rds_models(path = "./data/processed", html = TRUE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	cn.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	11186.547	18465.237	206225990
2	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 106.3)	cp.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	3800.096	3310.208	1179254390
3	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ Treatment + (1 \| quadrat)	20000	4	1	10000	5	10254.769	11806.319	224004760
4	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ n.treat + (1 \| quadrat)	20000	4	1	10000	17	8997.572	8327.074	956380015
5	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.litter ~ p.treat + (1 \| quadrat)	20000	4	1	10000	8	10806.876	12607.636	1236972411
6	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ n.treat + (1 \| quadrat)	20000	4	1	10000	7	8165.879	9671.253	1683655054
7	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ p.treat + (1 \| quadrat)	20000	4	1	10000	4	8190.692	11336.204	326631910
8	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	sil.k.wood ~ Treatment + (1 \| quadrat)	20000	4	1	10000	465	295.505	84.632	1147311351
9	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	litter.n.content.prop.initial ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	12453.012	17014.307	272268972
10	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.n ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	14603.282	16880.616	757288462
11	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	litter.p.content.prop.initial ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	8457.580	15360.879	2063397055
12	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.p ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	16204.903	17372.163	1797840914
13	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.n ~ prom.hoja.reman.sc * treat + (1 \| plot)	20000	4	1	10000	0	12014.609	18513.789	202030287
14	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 5.3)	np.mol.kg ~ treat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	11291.391	17616.689	226626333
15	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 2.5)	perc.p ~ prom.hoja.reman.sc * treat + (1 \| plot)	20000	4	1	10000	0	7081.653	10836.794	8818990
16	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	prop.litter.rem ~ trat * mo(days.fc) + (1 \| plot.f)	20000	4	1	10000	0	12245.978	13332.414	407718541
17	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.litter.rem ~ trat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	13040.090	18549.848	1631599382
18	b-normal(0, 10) Intercept-normal(0, 50) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 20)	prop.wood.rem ~ trat * days.sc + (1 \| plot.f)	20000	4	1	10000	0	11933.976	17029.142	1828689869

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brms_2.18.0       Rcpp_1.0.9        viridis_0.6.2     viridisLite_0.4.1
 [5] ggridges_0.5.4    posterior_1.3.1   ggrepel_0.9.1     cowplot_1.1.1    
 [9] tidybayes_3.0.1   ggplot2_3.4.0     readxl_1.3.1      knitr_1.42       
[13] kableExtra_1.3.4  klippy_0.0.0.9500

loaded via a namespace (and not attached):
  [1] backports_1.4.1      systemfonts_1.0.4    plyr_1.8.7          
  [4] igraph_1.3.5         splines_4.1.0        svUnit_1.0.6        
  [7] crosstalk_1.2.0      TH.data_1.1-0        rstantools_2.2.0    
 [10] inline_0.3.19        digest_0.6.31        htmltools_0.5.4     
 [13] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
 [16] remotes_2.4.2        RcppParallel_5.1.5   matrixStats_0.62.0  
 [19] xts_0.12.2           sandwich_3.0-1       svglite_2.1.0       
 [22] prettyunits_1.1.1    colorspace_2.0-3     rvest_1.0.3         
 [25] ggdist_3.2.0         xfun_0.36            dplyr_1.0.10        
 [28] callr_3.7.3          crayon_1.5.2         jsonlite_1.8.4      
 [31] lme4_1.1-27.1        survival_3.2-11      zoo_1.8-11          
 [34] glue_1.6.2           gtable_0.3.1         emmeans_1.8.1-1     
 [37] webshot_0.5.4        distributional_0.3.1 pkgbuild_1.4.0      
 [40] rstan_2.21.7         abind_1.4-5          scales_1.2.1        
 [43] mvtnorm_1.1-3        DBI_1.1.1            xaringanExtra_0.7.0 
 [46] miniUI_0.1.1.1       xtable_1.8-4         stats4_4.1.0        
 [49] StanHeaders_2.21.0-7 DT_0.26              htmlwidgets_1.5.4   
 [52] httr_1.4.4           threejs_0.3.3        arrayhelpers_1.1-0  
 [55] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1.9000      
 [58] farver_2.1.1         utf8_1.2.2           labeling_0.4.2      
 [61] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
 [64] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
 [67] tools_4.1.0          cli_3.6.0            generics_0.1.3      
 [70] evaluate_0.20        stringr_1.5.0        fastmap_1.1.0       
 [73] yaml_2.3.7           processx_3.8.0       purrr_1.0.0         
 [76] packrat_0.9.0        pbapply_1.6-0        nlme_3.1-152        
 [79] mime_0.12            projpred_2.0.2       formatR_1.11        
 [82] xml2_1.3.3           compiler_4.1.0       bayesplot_1.9.0     
 [85] shinythemes_1.2.0    rstudioapi_0.14      gamm4_0.2-6         
 [88] sketchy_1.0.2        tibble_3.1.8         stringi_1.7.12      
 [91] ps_1.7.2             Brobdingnag_1.2-9    lattice_0.20-44     
 [94] Matrix_1.5-1         nloptr_1.2.2.2       markdown_1.3        
 [97] shinyjs_2.1.0        tensorA_0.36.2       vctrs_0.5.2         
[100] pillar_1.8.1         lifecycle_1.0.3      bridgesampling_1.1-2
[103] estimability_1.4.1   httpuv_1.6.6         R6_2.5.1            
[106] promises_1.2.0.1     gridExtra_2.3        codetools_0.2-18    
[109] boot_1.3-28          colourpicker_1.2.0   MASS_7.3-54         
[112] gtools_3.9.3         assertthat_0.2.1     rprojroot_2.0.3     
[115] withr_2.5.0          shinystan_2.6.0      multcomp_1.4-17     
[118] mgcv_1.8-36          parallel_4.1.0       grid_4.1.0          
[121] minqa_1.2.4          tidyr_1.1.3          coda_0.19-4         
[124] rmarkdown_2.20       shiny_1.7.3          base64enc_0.1-3     
[127] dygraphs_1.1.1.6